Lectura y limpieza de Datos

Lectura de todos los archivos .csv de nuestro directorio link, previamente se hacen algunos ajustes a los datos originales, como quitar las comillas en los datos de 2021, dejar en una sola colmna la informaicón de 2022, entre otros. Y teniendo en cuenta la información dispuesta en cada uno de los archivos:

info_radares <- unique(rbind(as.matrix(flujo_vehicular_2019.csv %>% distinct(disp_nombre, lat, long)),
                             as.matrix(flujo_vehicular_2020.csv %>% distinct(disp_nombre, lat, long)),
                             as.matrix(flujo_vehicular_2022.csv %>% distinct(disp_nombre, lat, long))))


info_radares <- na.omit(data.frame(info_radares))

flujo_vehicular_2019 <- merge(x=flujo_vehicular_2019.csv,y=info_radares,by="disp_nombre")
names(flujo_vehicular_2019)[names(flujo_vehicular_2019) == "lat.y"] <- "lat"
names(flujo_vehicular_2019)[names(flujo_vehicular_2019) == "long.y"] <- "long"
flujo_vehicular_2019 <- flujo_vehicular_2019 %>% distinct(fecha, hora, autopista_nombre,
                                                          disp_nombre, disp_ubicacion,seccion_sentido,
                                                          lat,long, cantidad)

flujo_vehicular_2020 <- merge(x=flujo_vehicular_2020.csv,y=info_radares,by="disp_nombre")
names(flujo_vehicular_2020)[names(flujo_vehicular_2020) == "lat.y"] <- "lat"
names(flujo_vehicular_2020)[names(flujo_vehicular_2020) == "long.y"] <- "long"
flujo_vehicular_2020 <- flujo_vehicular_2020 %>% distinct(fecha, hora, autopista_nombre,
                                                          disp_nombre, disp_ubicacion,seccion_sentido,
                                                          lat,long, cantidad)

flujo_vehicular_2021 <- merge(x=flujo_vehicular_2021.csv,y=info_radares,by="disp_nombre")
names(flujo_vehicular_2021)[names(flujo_vehicular_2021) == "lat.y"] <- "lat"
names(flujo_vehicular_2021)[names(flujo_vehicular_2021) == "long.y"] <- "long"
flujo_vehicular_2021 <- flujo_vehicular_2021 %>% distinct(fecha, hora, autopista_nombre,
                                                          disp_nombre, disp_ubicacion,seccion_sentido,
                                                          lat,long, cantidad)

flujo_vehicular_2022 <- merge(x=flujo_vehicular_2022.csv,y=info_radares,by="disp_nombre")
names(flujo_vehicular_2022)[names(flujo_vehicular_2022) == "lat.y"] <- "lat"
names(flujo_vehicular_2022)[names(flujo_vehicular_2022) == "long.y"] <- "long"
flujo_vehicular_2022 <- flujo_vehicular_2022 %>% distinct(fecha, hora, autopista_nombre,
                                                          disp_nombre, disp_ubicacion,seccion_sentido,
                                                          lat,long, cantidad)


rm(flujo_vehicular_2019.csv,flujo_vehicular_2020.csv,flujo_vehicular_2021.csv,flujo_vehicular_2022.csv)
df <- unique(rbind(flujo_vehicular_2019,
                   flujo_vehicular_2020,
                   flujo_vehicular_2021,
                   flujo_vehicular_2022))

#write.table(df, file = "flujo_vehicular.txt", sep = "\t",
#            row.names = TRUE, col.names = NA)

df[, 'seccion_sentido'] <- as.factor(df[, 'seccion_sentido'])
df[, 'lat'] <- as.double(df[, 'lat'])
df[, 'long'] <- as.double(df[, 'long'])
df[, 'fecha'] <- as.Date((df[, 'fecha']), "%Y-%m-%d") 

summary(df)
##  disp_nombre            fecha                 hora       autopista_nombre  
##  Length:730833      Min.   :2019-01-01   Min.   : 0.00   Length:730833     
##  Class :character   1st Qu.:2019-10-31   1st Qu.: 5.00   Class :character  
##  Mode  :character   Median :2020-08-31   Median :11.00   Mode  :character  
##                     Mean   :2020-08-17   Mean   :11.47                     
##                     3rd Qu.:2021-06-16   3rd Qu.:17.00                     
##                     Max.   :2022-02-28   Max.   :23.00                     
##  disp_ubicacion  seccion_sentido    cantidad          lat        
##  Min.   :0.000    :  1045        Min.   :    0   Min.   :-34.68  
##  1st Qu.:1.050   A:362818        1st Qu.:  707   1st Qu.:-34.65  
##  Median :3.260   B:366970        Median : 1972   Median :-34.64  
##  Mean   :3.844                   Mean   : 2417   Mean   :-34.62  
##  3rd Qu.:7.200                   3rd Qu.: 3597   3rd Qu.:-34.55  
##  Max.   :9.900                   Max.   :19936   Max.   :-34.54  
##       long       
##  Min.   :-58.48  
##  1st Qu.:-58.47  
##  Median :-58.45  
##  Mean   :-58.43  
##  3rd Qu.:-58.38  
##  Max.   :-58.37

Después de consolidados los datos en una sola tabla, las fechas son transformadas en un solo formtao único para todos los registros, y se generan dos nuevas bases:

Los flujos vehiculares diarios se utilizan para algunas representaciones gráficas y/o mapas, los semanales se cargan a continuación para hacer algunos análisis adicionales y estimaciones.

my_data <- read_excel("flujo_vehicular_wcm.xlsx")
head(my_data, 5)
## # A tibble: 5 × 8
##   disp_nombre   WCM                 autopista_nombre disp_ubicacion
##   <chr>         <dttm>              <chr>                     <dbl>
## 1 RD100 Montiel 2019-07-22 00:00:00 AU Dellepiane               0.4
## 2 RD100 Montiel 2019-11-25 00:00:00 AU Dellepiane               0.4
## 3 RD100 Montiel 2019-11-04 00:00:00 AU Dellepiane               0.4
## 4 RD100 Montiel 2019-09-23 00:00:00 AU Dellepiane               0.4
## 5 RD100 Montiel 2019-09-16 00:00:00 AU Dellepiane               0.4
## # … with 4 more variables: seccion_sentido <chr>, lat <dbl>, long <dbl>,
## #   cantidad <dbl>

Ejercicio 1.

Genere una exploración de datos que permita la comprensión de los mismos (gráficos, medidas de tendencia central, mapas). Tenga en cuenta que muchas personas de negocio mirarán su informe por lo cual deberá facilitarles la comprensión del estado de situación.

Dado que no disponemos de muchas variables, algunas por las cuales podemos hacer apertura por categorías son, días de la semana, autopistas, y sección sentido, variables para las cuales harmos un primer acercamiento:

Se tienen en total cinco autopistas diferentes en la base de datos de flufo vehícular, sobre los cuales se puede concluir lo siguiente:

Analizando la información promedio por día de la semana, se observa

En cuanto al sentido de las autopistas, se tiene que:

Dashboard PBI

A continuación se presentan algunas imágenes tomadas del dashboard interactivo construido en la herramienta Power BI, donde se utilizaron mapas de ARCGis para tener una representación geoespacial de donde se encuentran los radares en la ciudad de Buenos Aires a lo largo de las autopistas analizadas anteriormente.

Mapa con informacíon promedio diaria de flujo vehicular para todos los radares

Ya es posible observar la ubicación de las autopistas y de los radares y concluir lo siguiente:

Lugones y Cantilo

Mapa con informacíon promedio diaria de flujo vehicular Cantilo & Lugones

En estas autopistas se tienen 8 de los 10 radares, donde se presenta mayor flujo vehícular; los radares Monroe, Esma, Udaondo, Dorrego y Ombues todos en sentido A sobre Lugones, son los radares donde se presenta mayor flujo vehicular en promedio por día.

Av9 julio sur

Mapa con informacíon promedio diaria de flujo vehicular Av9_julio_sur

El radar con mayor flujo vehicular de esta zona es Samperio A sobre la avenida 9 de Julio al sur oriente de la Ciudad

Ejercicio 2.

Desarrolle algún modelo que permita estimar el flujo vehicular para las próximas semanas. Justifique la elección del mismo e indique bajo qué medidas evaluaría su performance (no es necesario que tenga una alta precisión, lo importante es la justificación de por qué lo eligió).

Flujo vehicular en radar Monroe & Análisis Exploratorio de Datos (EDA)

Análisis de autocorrelación para examinar la dependencia en serie. Se utiliza para estimar qué valor en el pasado tiene una correlación con el valor actual. Proporciona la estimación p, d, q para los modelos ARIMA (Autoregresivo Integrado de Media Móvil).Análisis espectral para examinar el comportamiento cíclico. Se realiza para describir cómo la variación en una serie temporal puede ser explicada por componentes cíclicos. También se conoce como un análisis en dominio de la frecuencia. Usando esto, los componentes periódicos en un ambiente ruidoso se pueden separar.Estimación y descomposición de la tendencia. Se utiliza para el ajuste estacional. Busca construir, a partir de una serie temporal observada, una serie de componentes (que podrían usarse para reconstruir la serie original) donde cada una de ellas tiene una característica determinada. Antes de realizar cualquier EDA en los datos, debemos comprender los tres componentes de los datos de una serie temporal:

En R podemos averiguar estos componentes de la serie temporal con los siguientes comandos:

my_data$V1=as.character(my_data$cantidad)
my_data$V1=as.numeric(my_data$V1)
data <- filter(my_data, disp_nombre == "RD169 Monroe" & WCM >= '2019-01-01' & WCM < '2022-02-27')
data <- data[order(data$WCM),]

# frecuencia semanal de nuestros datos
RD169_Monroe <- ts(data$V1,start=c(2019,01,01),frequency = 52)
ajuste <- auto.arima(y = RD169_Monroe)
ts_decompose(RD169_Monroe)
par(mfrow=c(1,2))
acf(RD169_Monroe,main="ACF serie flujos", lag.max=36)
pacf(RD169_Monroe,main="PACF serie flujos", lag.max=36)

ACF y PACF. Correlogramas. Para saber cómo aplicar los modelos ARIMA tendremos que aprender a interpretar las gráficas de autocorrelación ACF y autocorrelación parcial PACF de una serie temporal. Estas gráficas nos ayundan a estimar los órdenes (recuerda…los 3 parámetros) (p,d,q) del modelo ARIMA.

Como regla, la PACF define el orden de AR(p) y la ACF el orden de MA(q). Otra cosa a tener encuentra para su interpretación es que hay que fijarse bien en los valores absolutos de los primeros datos, obviando las inversiones abajo-arriba que solo son indicadores de algún coeficiente negativo del modelo, pero no afectan al órden (p,q).

Si hay muchos valores de ACF altos(significativos), la serie es NO estacionaria, como en este caso, y habrá que agregar una diferenciación con d, y algo que podemos ver también en el test de Dickey-fuller, donde p-value > 0.05 y debemos aceptar la hipótesis nula de NO estacionaridad.

adf.test(RD169_Monroe)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  RD169_Monroe
## Dickey-Fuller = -2.1765, Lag order = 5, p-value = 0.5028
## alternative hypothesis: stationary
auto.arima(RD169_Monroe, trace=T)
## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2)            with drift         : 4043.409
##  ARIMA(0,1,0)            with drift         : 4070.493
##  ARIMA(1,1,0)            with drift         : 4061.47
##  ARIMA(0,1,1)            with drift         : 4043.534
##  ARIMA(0,1,0)                               : 4068.445
##  ARIMA(1,1,2)            with drift         : 4042.054
##  ARIMA(0,1,2)            with drift         : 4038.93
##  ARIMA(0,1,3)            with drift         : 4040.684
##  ARIMA(1,1,1)            with drift         : 4042.268
##  ARIMA(1,1,3)            with drift         : 4044.218
##  ARIMA(0,1,2)                               : 4036.877
##  ARIMA(0,1,1)                               : 4041.497
##  ARIMA(1,1,2)                               : 4039.951
##  ARIMA(0,1,3)                               : 4038.6
##  ARIMA(1,1,1)                               : 4040.184
##  ARIMA(1,1,3)                               : 4042.085
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,1,2)                               : 4061.423
## 
##  Best model: ARIMA(0,1,2)
## Series: RD169_Monroe 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.4384  -0.2153
## s.e.   0.0808   0.0800
## 
## sigma^2 = 2.284e+10:  log likelihood = -2027.63
## AIC=4061.26   AICc=4061.42   BIC=4070.33
#Se estiman posibles modelos
mod1<-arima(RD169_Monroe,order=c(0,1,2))
at_est <- residuals(mod1)
windows()
plot(at_est, main = "Residuos")

Ljung_Box <- NULL
x <- at_est
for(i in 1:18){
  Ljung_Box1 <- Box.test(x, lag = i, type="Ljung")$p.value
  Ljung_Box <- rbind(Ljung_Box, cbind(i,Ljung_Box1))
}
colnames(Ljung_Box) <- c("Rezago","p-valor");
cat(" \n \n")
cat("Prueba Ljung_Box Hipótesis nula independencia \n \n")
## Prueba Ljung_Box Hipótesis nula independencia 
## 
Ljung_Box
##       Rezago   p-valor
##  [1,]      1 0.8755701
##  [2,]      2 0.8744409
##  [3,]      3 0.9625602
##  [4,]      4 0.8875198
##  [5,]      5 0.9169196
##  [6,]      6 0.9257872
##  [7,]      7 0.9574069
##  [8,]      8 0.8843649
##  [9,]      9 0.7176691
## [10,]     10 0.7884379
## [11,]     11 0.7604969
## [12,]     12 0.8158561
## [13,]     13 0.8288200
## [14,]     14 0.8185509
## [15,]     15 0.8201282
## [16,]     16 0.7856981
## [17,]     17 0.8314543
## [18,]     18 0.8612310
plot(Ljung_Box, main="Prueba Ljung and Box \n Ho: Independencia",ylim=c(0,1))
abline(h=0.05,col="red")

### pruebas de hipótesis de normalidad
#ad.test(at_est) #test de normalidad (nortest)
shapiro.test(at_est) #test de normalidad (stats)
## 
##  Shapiro-Wilk normality test
## 
## data:  at_est
## W = 0.91647, p-value = 9.994e-08
win.graph(width=2.5,height=2.5,pointsize=4)
qqnorm(at_est); qqline(at_est,col=2)

win.graph(width=2.5,height=2.5,pointsize=8)
acf(at_est, main="ACF Residuos")

pacf(at_est)

Todos los valores de p para la prueba Ljung-Box Q están muy por encima de 0.05 , lo que nos indica que los datos no son dependientes.

Los valores son normales ya que descansan en una línea y no están por todos lados.

Como todos los gráficos apoyan el supuesto de que no hay un patrón en los residuales, podemos seguir adelante y calcular el pronóstico.

predicciones <- forecast(ajuste,h=5)
autoplot(predicciones)

predicciones[["mean"]]
## Time Series:
## Start = c(2021, 50) 
## End = c(2022, 2) 
## Frequency = 52 
## [1] 835393.0 833959.8 833959.8 833959.8 833959.8

ARIMA es una metodología de análisis y estimación para las series de tiempo, la cual es comunmente usada en estadística y econometría. Allí se tienen en cuenta aspectos muy importantes dichas series como lo son los procesos autoregresivos, las medias móviles, la estacionalidad y tendencia. En este caso de análisis, solamente se dispone de la variable temporal de flujo vehicular, por lo cual la metodología ARIMA es de correcta aplicación.

Pronóstico del flujo vehicular promedio semanal en Buenos Aires

En este caso haremos la estimación para teniendo en cuenta la totalidad de los radares

Serie de tiempo de la totalidad del flujo vehicular en la ciudad

library(dplyr)
df_total <- my_data %>% 
  dplyr::group_by(WCM) %>% 
  dplyr::summarise(mean_cantidad = mean(cantidad))

data <- filter(df_total, WCM >= '2019-01-01' &  WCM < '2022-02-27')
data <- filter(data,!is.na(mean_cantidad))
data <- data[order(df_total$WCM),]

total <- ts(data$mean_cantidad,start=c(2019,01,01),frequency = 52)
ajuste <- auto.arima(y = total)
ts_decompose(total)
predicciones <- forecast(ajuste,h=5)
autoplot(predicciones)
## Warning: Removed 2 row(s) containing missing values (geom_path).

predicciones[["mean"]]
## Time Series:
## Start = c(2022, 3) 
## End = c(2022, 7) 
## Frequency = 52 
## [1] 433702.9 425725.4 423553.3 426204.3 426543.8

En este caso particular tenemos comportamiento similar al revisado para la serie anterior que nos describía el comportamiento del radar Monroe en el tiempo, otro posible ajuste en nuestra metodología es cambiar la ventana de tiempo, y excluir inicios del 2020, datos que se vieron afectados por la pandemia.

data_2 <- filter(df_total, WCM >= '2020-05-01' & WCM < '2022-02-27')
data_2 <- filter(data_2,!is.na(mean_cantidad))
data_2 <- data_2[order(data_2$WCM),]

total <- ts(data_2$mean_cantidad,start=c(2020,05,01),frequency = 52)
ajuste <- auto.arima(y = total)

predicciones <- forecast(ajuste,h=5)
autoplot(predicciones)

predicciones[["mean"]]
## Time Series:
## Start = c(2021, 40) 
## End = c(2021, 44) 
## Frequency = 52 
## [1] 441022.5 444616.5 447603.1 450589.8 453576.4

Vemos por ejemplo, como la exclusión de un periodo de tiempo en particular puede cambiar nuestra estimación, y ahora vemos una ligera tendencia en aumento para las 5 semanas siguientes.

Ejercicio 3.

¿Cómo delegaría el monitoreo de este modelo en alguien del negocio? ¿Cuál cree que sería la mejor estrategia?

El pronóstico de la metodología viene acompañado de un intervalo de confianza, el cual se puede ajustar a la precisión que desee negocio o que se sugiera estadisticamente. Mensualmente se puede entregar un pronóstico de las semanas próximas con su respectivo intervalo de confianza, información que será util para realizar seguimiento de los datos reales. Dicha comparación permite generar alertas para la corecta toma de decisiones o acciones, e inclusive reentrenar la metodología.